2.1 Accessing the data in R

The pRolocdata package is a R Bioconductor experiment package that collects mass spectrometry-based spatial/organelle and protein-complex datasets from published experiments. Each data is stored in a container called a MSnSet instance (see the MSnbase for details) which allows users to work seamlessly with the R pRoloc and pRolocGUI software for spatial proteomics data analysis and visualisation.

The first step is to install the pRoloc and pRolocdata packages from Bioconductor

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("pRoloc", "pRolocdata")

## we also use the knitr and VennDiagram, ggplot2 packages in this 
## vignette so please  install if you do not have them 
BiocManager::install(c("knitr", "VennDiagram", "ggplot2", "dplyr", "kableExtra", 
                       "coda", "reshape2", "circlize", "ggalluvial", "tidyr", 
                       "Rmisc", "gridExtra", "gplots", "scico", "plotrix",
                       "RColorBrewer"))

Note: To access the widest selection on proteomics datasets we advise users to install the latest version pRolocdata package from Bioconductor. The THP data is in version >=1.29.1

Now we can load the packages by calling library function in R

library("pRoloc")
library("pRolocdata")

and subsequently we can now access all the datasets in pRolocdata. To list all available datasets we use the data function.

data(package = "pRolocdata")

We see there are 27 datasets (including the PSM and protein level data) associated with the manuscript [1], 20 from the hyperLOPIT data anlysis and 7 from the LPS timecourse.² For more information on these datasets, we can type

?thpLOPIT_lps_mulvey2021
?lpsTimecourse_mulvey2021

Figure 1. Screenshot of the documentation files in R for the THP LOPIT and LPS timecourse datasets

Each dataset can be loaded using the data function. For example, to load the final protein level LOPIT datasets for the unstimulated and 12 hour LPS stimulated datasets,

data("thpLOPIT_unstimulated_mulvey2021")
thpLOPIT_unstimulated_mulvey2021

data("thpLOPIT_lps_mulvey2021")
thpLOPIT_lps_mulvey2021

In the next sections we introduce the data, experimental design and brief overview of the experimental protocol before walking users through the downstream analysis which was used to generate the datasets available in pRolocdata and the manuscript Mulvey et al 2021.²

3.1 Summary

We assessed the global proteomic landscape of the THP-1 monocytic leukaemic cell line using a shotgun proteomics time course. Total (unfractionated) cell lysates were collected at 0, 2, 4, 6, 12, and 24 hour time-points following LPS stimulation and processed for quantitative, MS-based proteomics analysis. We reproducibly quantified 4,292 proteins across 3 biological replicate experiments, at all time-points. In total 311 proteins were found to be altered in abundance during the time course of LPS stimulation, with statistical significance (adjusted p-value < 0.01, log_2 fold-change > 0.6). The results are reported in full in Mulvey et al. 2021.²

3.2 PSM level data processing

The 3 PSM level datasets were imported into R and missing values were carefully assessed.

We can load the PSM and protein level data from the pRolocdata package by using the data function.

## Unstimulated data
data("psms_lpsTimecourse_rep1_mulvey2021")
data("psms_lpsTimecourse_rep2_mulvey2021")
data("psms_lpsTimecourse_rep3_mulvey2021")

psms_lpsTimecourse_rep1_mulvey2021

## MSnSet (storageMode: lockedEnvironment)
## assayData: 34640 features, 6 samples 
##   element names: exprs 
## protocolData: none
## phenoData
##   sampleNames: 0 hr rep 1 2 hr rep 1 ... 24 hr rep 1 (6 total)
##   varLabels: Tag Time Replicate
##   varMetadata: labelDescription
## featureData
##   featureNames: 1 2 ... 34640 (34640 total)
##   fvarLabels: Checked Confidence ... PSM.count (42 total)
##   fvarMetadata: labelDescription
## experimentData: use 'experimentData(object)'
## Annotation:  
## - - - Processing information - - -
## Subset [218958,6][34640,6] Tue Jan 12 14:42:27 2021 
##  MSnbase version: 2.14.2

dim(psms_lpsTimecourse_rep1_mulvey2021)

## [1] 34640     6

psms_lpsTimecourse_rep2_mulvey2021

## MSnSet (storageMode: lockedEnvironment)
## assayData: 49806 features, 6 samples 
##   element names: exprs 
## protocolData: none
## phenoData
##   sampleNames: 0 hr rep 2 2 hr rep 2 ... 24 hr rep 2 (6 total)
##   varLabels: Tag Time Replicate
##   varMetadata: labelDescription
## featureData
##   featureNames: 34641 34642 ... 84446 (49806 total)
##   fvarLabels: Checked Confidence ... PSM.count (42 total)
##   fvarMetadata: labelDescription
## experimentData: use 'experimentData(object)'
## Annotation:  
## - - - Processing information - - -
## Subset [218958,6][49806,6] Tue Jan 12 14:42:27 2021 
##  MSnbase version: 2.14.2

dim(psms_lpsTimecourse_rep2_mulvey2021)

## [1] 49806     6

psms_lpsTimecourse_rep3_mulvey2021

## MSnSet (storageMode: lockedEnvironment)
## assayData: 41853 features, 6 samples 
##   element names: exprs 
## protocolData: none
## phenoData
##   sampleNames: 0 hr rep 3 2 hr rep 3 ... 24 hr rep 3 (6 total)
##   varLabels: Tag Time Replicate
##   varMetadata: labelDescription
## featureData
##   featureNames: 84447 84448 ... 126299 (41853 total)
##   fvarLabels: Checked Confidence ... PSM.count (42 total)
##   fvarMetadata: labelDescription
## experimentData: use 'experimentData(object)'
## Annotation:  
## - - - Processing information - - -
## Subset [218958,6][41853,6] Tue Jan 12 14:42:27 2021 
##  MSnbase version: 2.14.2

dim(psms_lpsTimecourse_rep3_mulvey2021)

## [1] 41853     6

The summary above tells us that 34640, 49806 and 41853 PSMs were detected in replicates 1, 2 and 3, respectively.

In the below code chunk we put the replicates into a MSnSetList for ease of data processing.

psms <- MSnSetList(list(psms_lpsTimecourse_rep1_mulvey2021,
                        psms_lpsTimecourse_rep2_mulvey2021,
                        psms_lpsTimecourse_rep3_mulvey2021))

## Display number of PSMs for each protein with a MV
xx <- lapply(psms, function(z) which(apply(exprs(z), 1, anyNA)))
(numpsms <- sapply(seq(psms), function(y)
  {sapply(fData(psms[[y]])$Master.Protein.Accessions[xx[[y]]],
       function(z) length(which(fData(psms[[y]])$Master.Protein.Accessions 
                                == z)))}))

## [[1]]
##   P04179   P04179   P09914   O14879   P06307   Q8TCB0   P09601   P09601 
##        6        6       15       25        1        2       15       15 
##   O14879   P09913   P00973   P01584   P04075   Q9BYX4   O14879   P20591 
##       25       16        7        7       84       12       25       41 
## P29728-2   Q9BYX4   Q13501 
##       15       12       12 
## 
## [[2]]
##   O14879   P04179   P12956   P09913   O14879   P09601 P11388-4 
##       33       11       55       19       33       24       43 
## 
## [[3]]
##   P20591   P09914   Q13501   P20591   P20591   Q9BYX4   P20591   O14879 
##       68       21       14       68       68        9       68       30 
##   Q15833 Q96J88-2   P05362   P09913   Q9BYK8 Q96J88-2   O14879   P20592 
##       16        8       19       13       17        8       30       38

sapply(xx, length)

##  1  2  3 
## 19  7 16

## plot reps
plotMV <- function(r) {
  tochk <- unique(names(numpsms[[r]]))
  for (i in 1:length(tochk)) {
  mm <- exprs(psms[[r]])[which(fData(psms[[r]])$Master.Protein.Accessions == tochk[i]), , drop = FALSE]
  matplot(t(mm), type = "l", xlab = "", ylab = "m/z", lwd = 2, 
          main = tochk[i], col = "grey", axes = FALSE)
  axis(2)
  axis(1, at = 1:ncol(mm), labels = colnames(mm), las = 2)
  if (any(is.na(mm))) {
    ycoord <- sapply(1:ncol(mm), function(z) which(is.na(mm[, z])))
  .l <- sapply(ycoord, length)
  .xi <- which(.l > 0)
  xcoord <- unlist(lapply(1:length(.xi), function(z) rep(.xi[z], .l[.xi[z]])))
  ycoord <- unlist(ycoord)
  points(x = xcoord, y = ycoord, bg = "red", col = "black", pch = 24, cex = 1.5)}
  }
}

par(mfrow = c(4, 4))
plotMV(r = 1)

par(mfrow = c(2, 4))
plotMV(r = 2)

par(mfrow = c(3, 4))
plotMV(r = 3)

## Split the data by replicate
psms_rep1 <- psms[[1]]
psms_rep2 <- psms[[2]]
psms_rep3 <- psms[[3]]

## define psms which are missing at random - see ?impute
fData(psms_rep1)$randna <- FALSE
fData(psms_rep1)$randna[11237] <- TRUE      # P06307
fData(psms_rep1)$randna[26194] <- TRUE      # P04075

# remove missing this PSM for protein P12956 as it has 5 missing values out of the 6 timepoints
# remove PSM 48162 for protein P113388-4 as 42 other PSMs available and it has no clear trend
psms_rep2 <- psms_rep2[-c(31435, 48162), ]

## replicate 1 imputation
psms_rep1 <- impute(psms_rep1, "mixed", 
                    randna = fData(psms_rep1)$randna, 
                    mar = "knn",
                    mnar = "MinDet")

## replicate 2 imputation                
psms_rep2 <- impute(psms_rep2, "MinDet")

## replicate 3 imputation
psms_rep3 <- impute(psms_rep3, "MinDet")

## Re-create MSnSetList of the data
psms <- MSnSetList(list(rep1 = psms_rep1, rep2 = psms_rep2, rep3 = psms_rep3))

## Display number of PSMs for each protein with a MV
numpsms <- sapply(seq(psms), function(y)
  {sapply(fData(psms[[y]])$Master.Protein.Accessions[xx[[y]]],
       function(z) length(which(fData(psms[[y]])$Master.Protein.Accessions 
                                == z)))})

par(mfrow = c(4, 4))
plotMV(r = 1)

par(mfrow = c(2, 4))
plotMV(r = 2)

par(mfrow = c(3, 4))
plotMV(r = 3)

3.2.2 Data normalisation

The data was visually assessed by examining pairwise scatter plots and MAplots of the data and then normalised using vsn followed by diff.median.

## make all PSMs uppercase so that PSMs with different PTMs are not
## considered as different sequences
for (i in seq(psms)) 
  fData(psms@x[[i]])$Annotated.Sequence <- toupper(fData(psms@x[[i]])$Annotated.Sequence)


## combine the data and normalise
cmb <- MSnbase::combine(psms[[1]], psms[[2]])
cmb <- MSnbase::combine(cmb, psms[[3]])
dim(cmb)

## Let's have a look at the intensities
boxplot(log(exprs(cmb)))

## Try a few different normalisation approaches
cmb.vsn <- normalise(cmb, "vsn")
cmb.median <- normalise(cmb, "diff.median")
cmb.vsn.median <- normalise(cmb.vsn, "diff.median")
cmb.quantiles <- normalise(cmb, "quantiles")

## Plot the data
par(mfrow = c(2, 3), las = 2, oma = c(6, 1, 1, 1))
boxplot(log(exprs(cmb)), main = "method = untransformed", ylab = "log intensity")
boxplot(log(exprs(cmb.median)), main = "method = diff.median", ylab = "log intensity")
boxplot(exprs(cmb.vsn), main = "method = vsn", ylab = "log intensity")    # NOTE: VSN internally logs the data 
boxplot(log(exprs(cmb.vsn.median)), main = "method = vsn + diff.median", ylab = "log intensity")
boxplot(log(exprs(cmb.quantiles)), main = "method = quantiles", ylab = "log intensity")

## We decide to go with 'vsn' followed by 'diff.median'
cmb <- cmb.vsn.median

Figure X. Boxplots displaying the (log) intensities of PSMs after applying different normalisation methods available in MSnbase

psms1 <- filterNA(cmb[, 1:6])
psms2 <- filterNA(cmb[, 7:12])
psms3 <- filterNA(cmb[, 13:18])
psms <- MSnSetList(list(psms1, psms2, psms3))

## Aggregate to protein using median PSM per protein group
prots <- lapply(seq(psms), function(z)
               combineFeatures(psms[[z]], 
                               groupBy = fData(psms[[z]])$Master.Protein.Accession,
                               method = "median", verbose = FALSE))


## Do data fusion - create concatenated complete dataset
prots[[1]] <- updateFvarLabels(prots[[1]])
prots[[2]] <- updateFvarLabels(prots[[2]])
prots[[3]] <- updateFvarLabels(prots[[3]])
cmbprots <- MSnbase::combine(prots[[1]], prots[[2]])
cmbprots <- MSnbase::combine(cmbprots, prots[[3]])
cmbprots <- MSnbase::filterNA(cmbprots)
dim(cmbprots)

## [1] 4292   18

3.3 Protein level data processing

We load the individual protein level datasets from pRolocdata using the data function.

data("lpsTimecourse_rep1_mulvey2021")
data("lpsTimecourse_rep2_mulvey2021")
data("lpsTimecourse_rep3_mulvey2021")

lpsTimecourse_rep1_mulvey2021

## MSnSet (storageMode: lockedEnvironment)
## assayData: 5221 features, 6 samples 
##   element names: exprs 
## protocolData: none
## phenoData
##   sampleNames: X0.hr.rep.1 X2.hr.rep.1 ... X24.hr.rep.1 (6 total)
##   varLabels: Tag Time Replicate
##   varMetadata: labelDescription
## featureData
##   featureNames: A0AVT1 A0FGR8-2 ... Q9Y6Y8 (5221 total)
##   fvarLabels: Checked.rep1 Confidence.rep1 ... Peptides.count.rep1 (44
##     total)
##   fvarMetadata: labelDescription
## experimentData: use 'experimentData(object)'
## Annotation:  
## - - - Processing information - - -
##  MSnbase version: 2.14.2

dim(lpsTimecourse_rep1_mulvey2021)

## [1] 5221    6

lpsTimecourse_rep2_mulvey2021

## MSnSet (storageMode: lockedEnvironment)
## assayData: 5812 features, 6 samples 
##   element names: exprs 
## protocolData: none
## phenoData
##   sampleNames: X0.hr.rep.2 X2.hr.rep.2 ... X24.hr.rep.2 (6 total)
##   varLabels: Tag Time Replicate
##   varMetadata: labelDescription
## featureData
##   featureNames: A0A0B4J2F0 A0AV96 ... Q9Y6Y8 (5812 total)
##   fvarLabels: Checked.rep2 Confidence.rep2 ... Peptides.count.rep2 (43
##     total)
##   fvarMetadata: labelDescription
## experimentData: use 'experimentData(object)'
## Annotation:  
## - - - Processing information - - -
##  MSnbase version: 2.14.2

dim(lpsTimecourse_rep2_mulvey2021)

## [1] 5812    6

lpsTimecourse_rep3_mulvey2021

## MSnSet (storageMode: lockedEnvironment)
## assayData: 5210 features, 6 samples 
##   element names: exprs 
## protocolData: none
## phenoData
##   sampleNames: X0.hr.rep.3 X2.hr.rep.3 ... X24.hr.rep.3 (6 total)
##   varLabels: Tag Time Replicate
##   varMetadata: labelDescription
## featureData
##   featureNames: A0AV96 A0AVT1 ... Q9Y6Y8 (5210 total)
##   fvarLabels: Checked.rep3 Confidence.rep3 ... Peptides.count.rep3 (43
##     total)
##   fvarMetadata: labelDescription
## experimentData: use 'experimentData(object)'
## Annotation:  
## - - - Processing information - - -
##  MSnbase version: 2.14.2

dim(lpsTimecourse_rep3_mulvey2021)

## [1] 5210    6

We see we have 5221, 5812 and 5210 proteins for replicates 1, 2 and 3, respectively.

cmb <- MSnbase::combine(lpsTimecourse_rep1_mulvey2021, 
                        lpsTimecourse_rep2_mulvey2021)
cmb <- MSnbase::combine(cmb, lpsTimecourse_rep3_mulvey2021)
cmb <- filterNA(cmb)
dim(cmb)

## [1] 4292   18

source("r/ggVolcano.R")
## Add Gene Names to the data for adding Gene Name labels when we plot the results
cmb <- addGeneName(cmb, fcol = "Protein.Descriptions.rep1", col.name = "GN")

3.4 Protein expression analysis

par(mfrow = c(1, 2))
plot2D(t(cmb), fcol = "Tag", main = "PCA: TMT tag")
plot2D(t(cmb), fcol = "Replicate", main = "PCA: Replicates")

source("r/ggVolcano.R")

library("limma")
## 0hr and 24hr
i <- grep("X0.hr", sampleNames(cmb))
j <- grep("X24.hr", sampleNames(cmb))
eset <- cmb[, c(i, j)]
data <- exprs(eset)
rownames(data) <- fData(cmb)[["GN"]]
head(data)

##           X0.hr.rep.1 X0.hr.rep.2 X0.hr.rep.3 X24.hr.rep.1 X24.hr.rep.2
## UBA6         4.579279    5.242608    4.844744     4.871262     5.282477
## ESYT2        4.875324    5.136520    4.496599     4.397845     4.994216
## UHRF1BP1L    4.488927    3.614179    4.463849     4.480803     3.777839
## SHTN1        4.926983    4.752739    4.299910     5.145240     4.807961
## ILVBL        4.936206    4.395639    3.690728     4.480962     4.126668
## SH3PXD2B     5.335262    3.672522    2.516002     7.062948     5.140535
##           X24.hr.rep.3
## UBA6          5.058791
## ESYT2         4.119505
## UHRF1BP1L     3.843625
## SHTN1         4.651524
## ILVBL         3.573552
## SH3PXD2B      3.916628

## Create a model matrix
Time <- as.factor(pData(eset)$Time)
Replicate <- as.factor(pData(eset)$Replicate)
(design <- model.matrix(~Replicate+Time))

##   (Intercept) Replicate2 Replicate3 Time24
## 1           1          0          0      0
## 2           1          1          0      0
## 3           1          0          1      0
## 4           1          0          0      1
## 5           1          1          0      1
## 6           1          0          1      1
## attr(,"assign")
## [1] 0 1 1 2
## attr(,"contrasts")
## attr(,"contrasts")$Replicate
## [1] "contr.treatment"
## 
## attr(,"contrasts")$Time
## [1] "contr.treatment"

## Perform the fit
fit <- lmFit(data, design)
fit <- eBayes(fit)

## Get the results
res_t24 <- topTable(fit, coef = "Time24", 
                 adjust.method="BH", number = Inf)
res_t24 <- addPlottingInfo(res_t24)
msnset_t24 <- addLimmaInfo(cmb, res_t24)

## Look at p-value distribution
hist(res_t24[, "P.Value"], breaks = 50, 
     main = "Histogram of p-values for 0/24hr LPS \n(limma's paired moderated t-test)", 
     xlab = "p-value")

## save data as a TSV/CSV
write.exprs(msnset_t24, fDataCols = 131:138,  sep = "\t", 
            file = "csv/timecourse_t24.tsv")

library("ggplot2")
library("dplyr")
library("ggrepel")
ggVolcano(res_t24, "24 hours", N = 20, lfc = 0.6, p = 0.01, 
          addLegend = TRUE, text.size = 12, base.text.size = 12, 
          legend.text.size = 12, geom.point.size = 2, 
          label.text.size = 4)

## 0hr and 12hr
i <- grep("X0.hr", sampleNames(cmb))
j <- grep("X12.hr", sampleNames(cmb))
eset <- cmb[, c(i, j)]
data <- exprs(eset)
rownames(data) <- fData(cmb)[["GN"]]
head(data)

##           X0.hr.rep.1 X0.hr.rep.2 X0.hr.rep.3 X12.hr.rep.1 X12.hr.rep.2
## UBA6         4.579279    5.242608    4.844744     4.584263     5.190164
## ESYT2        4.875324    5.136520    4.496599     4.594140     5.180463
## UHRF1BP1L    4.488927    3.614179    4.463849     4.505748     3.814994
## SHTN1        4.926983    4.752739    4.299910     4.749325     4.655479
## ILVBL        4.936206    4.395639    3.690728     4.747260     4.037497
## SH3PXD2B     5.335262    3.672522    2.516002     6.088538     3.257607
##           X12.hr.rep.3
## UBA6          4.634475
## ESYT2         4.501894
## UHRF1BP1L     4.022813
## SHTN1         4.252117
## ILVBL         3.502920
## SH3PXD2B      2.956236

## Create a model matrix
Time <- as.factor(pData(eset)$Time)
Replicate <- as.factor(pData(eset)$Replicate)
design <- model.matrix(~Replicate+Time)

## Perform the fit
fit <- lmFit(data, design)
fit <- eBayes(fit)

## Get the results
res_t12 <- topTable(fit, coef = "Time12", 
                 adjust.method="BH", number = Inf)
res_t12 <- addPlottingInfo(res_t12)
msnset_t12 <- addLimmaInfo(cmb, res_t12)

## Look at p-value distribution
hist(res_t12[, "P.Value"], breaks = 50, 
     main = "Histogram of p-values for 0/12hr LPS \n(limma's paired moderated t-test)", 
     xlab = "p-value")

## save data as a TSV/CSV
write.exprs(msnset_t12, fDataCols = 131:138,  sep = "\t", 
            file = "csv/timecourse_t12.tsv")

## plot volcanos
ggVolcano(res_t12, "12 hours", N = 20, lfc = 0.6, 
          p = 0.01, addLegend = FALSE,
          text.size = 12, base.text.size = 12, 
          legend.text.size = 12, geom.point.size = 2, 
          label.text.size = 4)

## Warning: Removed 12 rows containing missing values (geom_label_repel).

## 0hr and 6r
i <- grep("X0.hr", sampleNames(cmb))
j <- grep("X6.hr", sampleNames(cmb))
eset <- cmb[, c(i, j)]
data <- exprs(eset)
rownames(data) <- fData(cmb)[["GN"]]
head(data)

##           X0.hr.rep.1 X0.hr.rep.2 X0.hr.rep.3 X6.hr.rep.1 X6.hr.rep.2
## UBA6         4.579279    5.242608    4.844744    4.478218    5.096956
## ESYT2        4.875324    5.136520    4.496599    4.913246    5.250985
## UHRF1BP1L    4.488927    3.614179    4.463849    4.546953    4.030109
## SHTN1        4.926983    4.752739    4.299910    4.968413    4.618690
## ILVBL        4.936206    4.395639    3.690728    4.783785    4.140193
## SH3PXD2B     5.335262    3.672522    2.516002    5.255045    3.250098
##           X6.hr.rep.3
## UBA6         4.839907
## ESYT2        4.330492
## UHRF1BP1L    4.313316
## SHTN1        4.708527
## ILVBL        4.019145
## SH3PXD2B     2.282174

## Create a model matrix
Time <- as.factor(pData(eset)$Time)
Replicate <- as.factor(pData(eset)$Replicate)
design <- model.matrix(~Replicate+Time)

## Perform the fit
fit <- lmFit(data, design)
fit <- eBayes(fit)

## Get the results
res_t6 <- topTable(fit, coef = "Time6", 
                 adjust.method="BH", number = Inf)
res_t6 <- addPlottingInfo(res_t6)
msnset_t6 <- addLimmaInfo(cmb, res_t6)

## Look at p-value distribution
hist(res_t6[, "P.Value"], breaks = 50, 
     main = "Histogram of p-values for 0/6hr LPS \n(limma's paired moderated t-test)", 
     xlab = "p-value")

## save data as a TSV/CSV
write.exprs(msnset_t6, fDataCols = 131:138,  sep = "\t", 
            file = "csv/timecourse_t6.tsv")

## plot volcanos
ggVolcano(res_t6, "6 hours", N = 20, lfc = 0.6, 
          p = 0.01, addLegend = FALSE,
          text.size = 12, base.text.size = 12, 
          legend.text.size = 12, geom.point.size = 2, 
          label.text.size = 4)

## 0hr and 4hr
i <- grep("X0.hr", sampleNames(cmb))
j <- grep("X4.hr", sampleNames(cmb))
eset <- cmb[, c(i, j)]
data <- exprs(eset)
rownames(data) <- fData(cmb)[["GN"]]
head(data)

##           X0.hr.rep.1 X0.hr.rep.2 X0.hr.rep.3 X4.hr.rep.1 X4.hr.rep.2
## UBA6         4.579279    5.242608    4.844744    4.674778    5.360775
## ESYT2        4.875324    5.136520    4.496599    4.993125    5.268153
## UHRF1BP1L    4.488927    3.614179    4.463849    4.604388    3.659608
## SHTN1        4.926983    4.752739    4.299910    5.127168    4.831314
## ILVBL        4.936206    4.395639    3.690728    4.810118    4.155659
## SH3PXD2B     5.335262    3.672522    2.516002    5.334079    3.722221
##           X4.hr.rep.3
## UBA6         4.844971
## ESYT2        4.744458
## UHRF1BP1L    4.296213
## SHTN1        4.886162
## ILVBL        3.649692
## SH3PXD2B     2.313917

## Create a model matrix
Time <- as.factor(pData(eset)$Time)
Replicate <- as.factor(pData(eset)$Replicate)
design <- model.matrix(~Replicate+Time)

## Perform the fit
fit <- lmFit(data, design)
fit <- eBayes(fit)

## Get the results
res_t4 <- topTable(fit, coef = "Time4", 
                 adjust.method="BH", number = Inf)
res_t4 <- addPlottingInfo(res_t4)
msnset_t4 <- addLimmaInfo(cmb, res_t4)

## Look at p-value distribution
hist(res_t4[, "P.Value"], breaks = 50, 
     main = "Histogram of p-values for 0/4hr LPS \n(limma's paired moderated t-test)", 
     xlab = "p-value")

## save data as a TSV/CSV
write.exprs(msnset_t4, fDataCols = c(131:138),  sep = "\t", 
            file = "csv/timecourse_t4.tsv")

## plot volcanos
ggVolcano(res_t4, "4 hours", N = 20, lfc = 0.6, 
          p = 0.01, addLegend = FALSE,
          text.size = 12, base.text.size = 12, 
          legend.text.size = 12, geom.point.size = 2, 
          label.text.size = 4)

## 0hr and 2hr
i <- grep("X0.hr", sampleNames(cmb))
j <- grep("X2.hr", sampleNames(cmb))
eset <- cmb[, c(i, j)]
data <- exprs(eset)
rownames(data) <- fData(cmb)[["GN"]]
head(data)

##           X0.hr.rep.1 X0.hr.rep.2 X0.hr.rep.3 X2.hr.rep.1 X2.hr.rep.2
## UBA6         4.579279    5.242608    4.844744    4.359477    5.190733
## ESYT2        4.875324    5.136520    4.496599    4.674483    5.146304
## UHRF1BP1L    4.488927    3.614179    4.463849    4.478383    3.873257
## SHTN1        4.926983    4.752739    4.299910    4.891723    4.652700
## ILVBL        4.936206    4.395639    3.690728    4.543970    4.174060
## SH3PXD2B     5.335262    3.672522    2.516002    5.101335    4.315416
##           X2.hr.rep.3
## UBA6         4.691351
## ESYT2        4.753582
## UHRF1BP1L    4.712662
## SHTN1        4.917707
## ILVBL        3.817152
## SH3PXD2B     2.197979

## Create a model matrix
Time <- as.factor(pData(eset)$Time)
Replicate <- as.factor(pData(eset)$Replicate)
design <- model.matrix(~Replicate+Time)

## Perform the fit
fit <- lmFit(data, design)
fit <- eBayes(fit)

## Get the results
res_t2 <- topTable(fit, coef = "Time2", 
                 adjust.method="BH", number = Inf)
res_t2 <- addPlottingInfo(res_t2)
msnset_t2 <- addLimmaInfo(cmb, res_t2)

## Look at p-value distribution
hist(res_t2[, "P.Value"], breaks = 50, 
     main = "Histogram of p-values for 0/2hr LPS \n(limma's paired moderated t-test)", 
     xlab = "p-value")

## save data as a TSV/CSV
write.exprs(msnset_t2, fDataCols = 131:138,  sep = "\t", 
            file = "csv/timecourse_t2.tsv")

## plot volcanos
ggVolcano(res_t2, "2 hours", N = 20, lfc = 0.6, 
          p = 0.05, addLegend = FALSE,
          text.size = 12, base.text.size = 12, 
          legend.text.size = 12, geom.point.size = 2, 
          label.text.size = 4)

source("r/plotLineGraph.R")
library("tidyr")
library("Rmisc")
library("gridExtra")

## read in elisa data
elisa_il1b_wide <- read.csv("csv/elisa_il1b.csv")
elise_cxcl_wide <- read.csv("csv/elisa_cscl10.csv")

## convert to long tidy format for ggplot
elisa_il1b_long <- gather(elisa_il1b_wide, replicate, measurement, X1:X6, factor_key = TRUE)
elisa_cxcl10_long <- gather(elise_cxcl_wide, replicate, measurement, X1:X6, factor_key = TRUE)

p <- 
  plotLineGraph(elisa_il1b_wide, "median", linecol = "#1B83E9") +
  ggtitle("Elisa data: IL-1B") + 
  theme(plot.title = element_text(hjust = 0.5))

q <- plotLineGraph(elise_cxcl_wide, "median", linecol = "#1B83E9") +
  ggtitle("Elisa data: CXCL10") + 
  theme(plot.title = element_text(hjust = 0.5))

## get quantitation data
IL1B <- "P01584"
CXCL10 <- "P02778"

## Load protein data
x1 <- lpsTimecourse_rep1_mulvey2021
x2 <- lpsTimecourse_rep2_mulvey2021
x3 <- lpsTimecourse_rep3_mulvey2021

## convert to long tidy format for ggplot
lopit_il1b_wide <- t(rbind(exprs(x1[IL1B,]), exprs(x2[IL1B,]),exprs(x3[IL1B,])))
lopit_cxcl10_wide <- t(rbind(exprs(x1[CXCL10,]), exprs(x2[CXCL10,]),exprs(x3[CXCL10,])))
lopit_il1b_wide <- data.frame(cbind(Time = c(0, 2, 4, 6, 12, 24), lopit_il1b_wide))
lopit_cxcl10_wide <- data.frame(cbind(Time = c(0, 2, 4, 6, 12, 24), lopit_cxcl10_wide))
colnames(lopit_cxcl10_wide) <- colnames(lopit_il1b_wide) <- c("Time", "X1", "X2", "X6")
lopit_il1b_long <- gather(lopit_il1b_wide, replicate, measurement, X1:X6, factor_key = TRUE)
lopit_cxcl10_long <- gather(lopit_cxcl10_wide, replicate, measurement, X1:X6, factor_key = TRUE)

## plot quantitation profiles
r = plotLineGraph(lopit_il1b_wide, 
                  linecol = "#C50000", 
                  .ylab = "Abundance\n",
                  themeSize = 14) +
  ylim(c(1.5, 7)) +
  ggtitle("Shotgun: IL-1B") + 
  theme(plot.title = element_text(hjust = 0.5))
s = plotLineGraph(lopit_cxcl10_wide, 
                  linecol = "#C50000", 
                  .ylab = "Abundance\n",
                  themeSize = 14) +
  ylim(c(1.5, 7)) +
  ggtitle("Shotgun: CXCL10") + 
  theme(plot.title = element_text(hjust = 0.5))

grid.arrange(p, q, r, s, ncol = 2, nrow = 2)

library("gplots")
library("plotrix")
library("scico")

## restructure the gather the data for heatmap
dat <- cbind(exprs(msnset_t12), fData(msnset_t12))
rownames(dat) <- NULL
dat <- dat %>% 
  filter(Significance != "Not significant") %>% 
  tibble::column_to_rownames("GN") %>% 
  select(1:18) %>% as.matrix()
calcMeds <- function(df, .grepName) {
  ind <- grep(.grepName, colnames(df))
  .df <- df[, ind]
  .mean <- rowMedians(.df)
  return(.mean)
}
t0 <- calcMeds(dat, "X0.hr")
t2 <- calcMeds(dat, "X2.hr")
t4 <- calcMeds(dat, "X4.hr")
t6 <- calcMeds(dat, "X6.hr")
t12 <- calcMeds(dat, "X12.hr")
t24 <- calcMeds(dat, "X24.hr")
dat_meds <- cbind(t0, t2, t4, t6, t12, t24)
rownames(dat_meds) <- rownames(dat)
colnames(dat_meds) <- c(0, 2, 4, 6, 12, 24)
head(dat_meds)

##              0        2        4        6       12       24
## IFIT3 2.730084 3.055495 4.169061 5.354609 6.142207 6.638291
## CETN3 5.825406 5.125166 5.315861 5.242166 4.666942 4.932849
## ZEB2  4.142148 4.718323 4.614038 4.773679 5.162687 5.100232
## ABCA1 3.824979 3.926887 4.264135 5.227943 5.590898 5.126276
## NADK  3.492452 3.460054 3.345293 3.647640 4.756142 5.302163
## DDX58 3.991000 3.916902 3.878920 4.849751 5.499738 5.950881

## set colour palette - use the scico package
my_palette <- scico(72, palette = "berlin")
col_breaks = c(seq(-1,1,length=35), 
               seq(0,0.8,length=11),  
               seq(1,0,length=35)) 

## Plot heatmap
heatmap.2(dat_meds, density.info = "none", trace = "none", 
          col = my_palette, Colv = NA, dendrogram = "row", 
          key=TRUE, cexCol = 2, cexRow = 1.3, margin = c(6,8), 
          lhei=c(.8,3.5), lwid=c(0.2,0.8), offsetRow=0.3)

3.5 Bayesian temporal clustering

Bayesian correlated clustering was conducted on the shotgun proteomic time-course series to capture clusters of co-regulated proteins and temporal patterns of protein expression in response to LPS. Replicate experiments were combined using multiple dataset integration MDI⁵ where the number of clusters were automatically inferred.

Bayesian inference was performed using Markov-chain Monte Carlo (MCMC) as implemented in the MDI-GPU software.⁵ Clusters were extracted only for proteins that were consistently allocated to the same clusterings (sampled from the posterior distribution) across replicates (the results can be found in Supplementary Table 3 of Mulvey et al. and also in the csv directory of this Git repo).

In the below code chunk we load the results from temporal clustering and plot 12 clusters of interest as discussed in Mulvey et al. (Figure 2 of²)

library("RColorBrewer")
source("r/plotRibbons.R")

## Read in MDI clustering results
mdi_rep1 <- read.csv("csv/mdi_clustering_lpsrep1.csv")
mdi_rep2 <- read.csv("csv/mdi_clustering_lpsrep2.csv")
mdi_rep3 <- read.csv("csv/mdi_clustering_lpsrep3.csv")
mdi <- read.csv("csv/mdi_clusters_output.csv")

# get median normalised intensity
dat_mat <- list(data.matrix(mdi_rep1[, -1]), 
                data.matrix(mdi_rep2[, -1]), 
                data.matrix(mdi_rep3[, -1]))
av_mat <- apply(simplify2array(dat_mat), c(1,2), median)
rownames(av_mat) <- mdi_rep1[,1]
colnames(av_mat) <- c(0, 2, 4, 6, 12, 24)
head(av_mat)

##                 0          2          4           6        12       24
## P02795 -1.0406225 -0.8355461 -0.6604120  0.13893606 1.0759237 1.258747
## P04179 -0.6899927 -0.6821170 -0.5093226 -0.30338106 0.1302568 1.952588
## O14879 -0.8874299 -0.8758429 -0.5680529 -0.01642289 0.6585623 1.661210
## Q9BYX4 -0.8292398 -0.8356126 -0.6969017 -0.08215123 0.9237005 1.515785
## P00973 -0.8814504 -0.8015016 -0.6501797 -0.11263572 0.8438697 1.578594
## P05362 -0.8188745 -0.7149888 -0.4673084 -0.22529338 0.5581240 1.756882

## plot clusters of interest
ids <- c(1:6, 9:11, 13, 14, 17)
mdi_profs <- sapply(ids, function(z) 
  av_mat[mdi$Protein.Accession[which(mdi$Cluster.Number == z)], ])
names(mdi_profs) <- paste0("cluster", ids)

## pick a good palette
ribbonCol <- c(brewer.pal(n = 9, "Oranges")[-c(1:3)],
           brewer.pal(n = 9, "GnBu")[-c(1:3)])

par(mfrow = c(3,2), mar = c(6, 5, 3, 3))
for (i in 1:6) {
  ribbonPlot(mdi_profs[[i]], col = ribbonCol[i], c(0.05, .95), 
             main = paste("Cluster", ids[i]), cex.main = 2, size = 2)
}

par(mfrow = c(3,2), mar = c(6, 5, 3, 3))
for (i in 7:12) {
  ribbonPlot(mdi_profs[[i]], col = ribbonCol[i], c(0.05, .95), 
             main = paste("Cluster", ids[i]), cex.main = 2, size = 2)
}

3.6 Gene Ontology enrichment

Temporal clusters were functionally annotated by using the Gene Ontology (GO) [xxx]. For each cluster output from the Bayesian temporal clustering, in turn, an enrichment test for biological process (BP), cellular component (CC) and molecular function (MF) annotations, was carried out using the Database for Annotation, Visualisation and Integrated Discovery (DAVID v6.8; https://david.ncifcrf.gov/) where a cutoff of < 0.05 was applied according to the Benjamini-Hochberg procedure [xxx]. The GO annotation enrichment results are found in Supplementary Table 4 of Mulvey et al 2021.

In the below code chunk some of the results from the GO enrichment (Figure 3 of Mulvey et al. 2021). GOCC, GOMF and GOBP annotation term enrichment for the the clusters are shown below where the -log10(adjusted p-value) is shown along the x-axis and most highly enriched term along the y-axis. The numbers within each barplot refer to the number of proteins associated with the term in each cluster.

source("r/GObarplot.R")
GObarplots(filename = "csv/david_data_gobp.csv", cols = ribbonCol, N = 30) +
  ggtitle("GO Biological Process Terms")

GObarplots(filename = "csv/david_data_gocc.csv", cols = ribbonCol) +
  ggtitle("GO Cellular Component Terms")

GObarplots(filename = "csv/david_data_gomf.csv", cols = ribbonCol, N = 15) +
  ggtitle("GO Molecular Function Terms")

The timecourse of LPS provides an insight into the protein expression of THP-1 proteome during the initial phases of a pro-inflammatory response. In order the gain a complete picture of protein relocalisation events we performed a series of hyperLOPIT experiments and 0 and 12 hour LPS time-points to give a deeper spaial understanding of the dynamic changes that occur in the proteome following stimulation.

4.1 The hyperLOPIT workflow

Using the hyperLOPIT proteomics platform⁶ we obtained spatial maps that capture the steady-state distribution of thousands of proteins in the THP-1 human leukaemia cell line. Experiments were conducted in triplicate and the samples were analysed and collected when cells were (1) unstimulated and then (2) at 12 hours following stimulation with LPS.

The hyperLOPIT method combines biochemical cell fractionation with multiplexed high-resolution mass-spectrometry,⁶ the full protocol and experimental design can be found in.² Briefly, TMT labelling was conducted as described in⁷ and 20 fractions including the cytosol-enriched fraction were labelled per gradient, by combining two TMT 10-plex experiments in an interleaved labelling design to capture as much subcellular diversity as possible. The experimental design for each dataset is stored in the pData slot of each dataset.

To access the full experimental design use the pData function

## TMT experiment data for the unstimulated hyperLOPIT data
pData(thpLOPIT_unstimulated_mulvey2021)

##                               Tag    Treatment Replicate Set Fraction
## unstim_rep1_set1_126_cyto     126 Unstimulated         1   1     cyto
## unstim_rep2_set1_126_cyto     126 Unstimulated         2   1     cyto
## unstim_rep3_set1_126_cyto     126 Unstimulated         3   1     cyto
## unstim_rep1_set1_127N_F1.4   127N Unstimulated         1   1     F1.4
## unstim_rep2_set1_127N_F1.6   127N Unstimulated         2   1     F1.6
## unstim_rep3_set1_127N_F1.6   127N Unstimulated         3   1     F1.6
## unstim_rep1_set2_126_F5.6     126 Unstimulated         1   2     F5.6
## unstim_rep1_set1_127C_F7     127C Unstimulated         1   1       F7
## unstim_rep2_set2_127N_F7.8   127N Unstimulated         2   2     F7.8
## unstim_rep3_set2_126_F7.8     126 Unstimulated         3   2     F7.8
## unstim_rep1_set2_127N_F8     127N Unstimulated         1   2       F8
## unstim_rep1_set1_128N_F9     128N Unstimulated         1   1       F9
## unstim_rep3_set1_127C_F9     127C Unstimulated         3   1       F9
## unstim_rep2_set1_127C_F9.10  127C Unstimulated         2   1    F9.10
## unstim_rep1_set2_127C_F10    127C Unstimulated         1   2      F10
## unstim_rep3_set2_127N_F10    127N Unstimulated         3   2      F10
## unstim_rep1_set1_128C_F11    128C Unstimulated         1   1      F11
## unstim_rep3_set1_128N_F11    128N Unstimulated         3   1      F11
## unstim_rep2_set2_127C_F11.12 127C Unstimulated         2   2   F11.12
## unstim_rep1_set2_128N_F12    128N Unstimulated         1   2      F12
## unstim_rep2_set1_128N_F12    128N Unstimulated         2   1      F12
## unstim_rep3_set2_127C_F12    127C Unstimulated         3   2      F12
## unstim_rep1_set1_129N_F13    129N Unstimulated         1   1      F13
## unsti_rep2_set2_128N_F13     128N Unstimulated         2   2      F13
## unstim_rep3_set1_128C_F13    128C Unstimulated         3   1      F13
## unstim_rep1_set2_128C_F14    128C Unstimulated         1   2      F14
## unstim_rep2_set1_128C_F14    128C Unstimulated         2   1      F14
## unstim_rep3_set2_128N_F14    128N Unstimulated         3   2      F14
## unstim_rep1_set1_129C_F15    129C Unstimulated         1   1      F15
## unstim_rep2_set2_128C_F15    128C Unstimulated         2   2      F15
## unstim_rep3_set1_129N_F15    129N Unstimulated         3   1      F15
## unstim_rep1_set2_129N_F16    129N Unstimulated         1   2      F16
## unstim_rep2_set1_129N_F16    129N Unstimulated         2   1      F16
## unstim_rep3_set2_128C_F16    128C Unstimulated         3   2      F16
## unstim_rep1_set1_130N_F17    130N Unstimulated         1   1      F17
## unstim_rep2_set2_129N_F17    129N Unstimulated         2   2      F17
## unstim_rep3_set1_129C_F17    129C Unstimulated         3   1      F17
## unstim_rep1_set2_129C_F18    129C Unstimulated         1   2      F18
## unstim_rep2_set1_129C_F18    129C Unstimulated         2   1      F18
## unstim_rep3_set2_129N_F18    129N Unstimulated         3   2      F18
## unstim_rep1_set1_130C_F19    130C Unstimulated         1   1      F19
## unstim_rep2_set2_129C_F19    129C Unstimulated         2   2      F19
## unstim_rep3_set1_130N_F19    130N Unstimulated         3   1      F19
## unstim_rep1_set2_130N_F20    130N Unstimulated         1   2      F20
## unstim_rep2_set1_130N_F20    130N Unstimulated         2   1      F20
## unstim_rep3_set2_129C_F20    129C Unstimulated         3   2      F20
## unstim_rep1_set1_131_F21      131 Unstimulated         1   1      F21
## unstim_rep2_set2_130N_F21    130N Unstimulated         2   2      F21
## unstim_rep3_set1_130C_F21    130C Unstimulated         3   1      F21
## unstim_rep1_set2_130C_F22    130C Unstimulated         1   2      F22
## unstim_rep2_set1_130C_F22    130C Unstimulated         2   1      F22
## unstim_rep3_set2_130N_F22    130N Unstimulated         3   2      F22
## unstim_rep2_set2_130C_F23    130C Unstimulated         2   2      F23
## unstim_rep3_set1_131_F23      131 Unstimulated         3   1      F23
## unstim_rep1_set2_131_F24      131 Unstimulated         1   2      F24
## unstim_rep2_set1_131_F24      131 Unstimulated         2   1      F24
## unstim_rep3_set2_130C_F24    130C Unstimulated         3   2      F24
## unstim_rep2_set2_131_F25      131 Unstimulated         2   2      F25
## unstim_rep3_set2_131_F26      131 Unstimulated         3   2      F26

## TMT labelling scheme for the LPS hyperLOPIT data
pData(thpLOPIT_lps_mulvey2021)

##                           Tag Treatment Replicate Set Fraction
## lps_rep1_set1_126_cyto    126       LPS         1   1     cyto
## lps_rep2_set1_126_cyto    126       LPS         2   1     cyto
## lps_rep3_set1_126_cyto    126       LPS         3   1     cyto
## lps_rep1_set1_127N_F1.4  127N       LPS         1   1     F1.4
## lps_rep2_set1_127N_F1.6  127N       LPS         2   1     F1.6
## lps_rep3_set1_127N_F1.6  127N       LPS         3   1     F1.6
## lps_rep1_set2_126_F5.7    126       LPS         1   2     F5.7
## lps_rep2_set2_126_F7.8    126       LPS         2   2     F7.8
## lps_rep3_set2_126_F7.8    126       LPS         3   2     F7.8
## lps_rep1_set1_127C_F8    127C       LPS         1   1       F8
## lps_rep1_set2_127N_F9    127N       LPS         1   2       F9
## lps_rep3_set1_127C_F9    127C       LPS         3   1       F9
## lps_rep2_set1_127C_F9.10 127C       LPS         2   1    F9.10
## lps_rep1_set1_128N_F10   128N       LPS         1   1      F10
## lps_rep3_set2_127N_F10   127N       LPS         3   2      F10
## lps_rep1_set2_127C_F11   127C       LPS         1   2      F11
## lps_rep2_set2_127N_F11   127N       LPS         2   2      F11
## lps_rep3_set1_128N_F11   128N       LPS         3   1      F11
## lps_rep1_set1_128C_F12   128C       LPS         1   1      F12
## lps_rep2_set1_128N_F12   128N       LPS         2   1      F12
## lps_rep3_set2_127C_F12   127C       LPS         3   2      F12
## lps_rep1_set2_128N_F13   128N       LPS         1   2      F13
## lps_rep2_set2_127C_F13   127C       LPS         2   2      F13
## lps_rep3_set1_128C_F13   128C       LPS         3   1      F13
## lps_rep1_set1_129N_F14   129N       LPS         1   1      F14
## lps_rep2_set1_128C_F14   128C       LPS         2   1      F14
## lps_rep3_set2_128N_F14   128N       LPS         3   2      F14
## lps_rep1_set2_128C_F15   128C       LPS         1   2      F15
## lps_rep2_set2_128N_F15   128N       LPS         2   2      F15
## lps_rep3_set1_129N_F15   129N       LPS         3   1      F15
## lps_rep1_set1_129C_F16   129C       LPS         1   1      F16
## lps_rep2_set1_129N_F16   129N       LPS         2   1      F16
## lps_rep3_set2_128C_F16   128C       LPS         3   2      F16
## lps_rep1_set2_129N_F17   129N       LPS         1   2      F17
## lps_rep2_set2_128C_F17   128C       LPS         2   2      F17
## lps_rep3_set1_129C_F17   129C       LPS         3   1      F17
## lps_rep1_set1_130N_F18   130N       LPS         1   1      F18
## lps_rep2_set1_129C_F18   129C       LPS         2   1      F18
## lps_rep3_set2_129N_F18   129N       LPS         3   2      F18
## lps_rep1_set2_129C_F19   129C       LPS         1   2      F19
## lps_rep2_set2_129N_F19   129N       LPS         2   2      F19
## lps_rep3_set1_130N_F19   130N       LPS         3   1      F19
## lps_rep1_set1_130C_F20   130C       LPS         1   1      F20
## lps_rep2_set1_130N_F20   130N       LPS         2   1      F20
## lps_rep3_set2_129C_F20   129C       LPS         3   2      F20
## lps_rep1_set2_130N_F21   130N       LPS         1   2      F21
## lps_rep2_set2_129C_F21   129C       LPS         2   2      F21
## lps_rep3_set1_130C_F21   130C       LPS         3   1      F21
## lps_rep1_set1_131_F22     131       LPS         1   1      F22
## lps_rep2_set1_130C_F22   130C       LPS         2   1      F22
## lps_rep3_set2_130N_F22   130N       LPS         3   2      F22
## lps_rep1_set2_130C_F23   130C       LPS         1   2      F23
## lps_rep2_set2_130N_F23   130N       LPS         2   2      F23
## lps_rep3_set1_131_F23     131       LPS         3   1      F23
## lps_rep2_set1_131_F24     131       LPS         2   1      F24
## lps_rep3_set2_130C_F24   130C       LPS         3   2      F24
## lps_rep1_set2_131_F25     131       LPS         1   2      F25
## lps_rep2_set2_130C_F25   130C       LPS         2   2      F25
## lps_rep3_set2_131_F26     131       LPS         3   2      F26

As described in⁷ the data was processed with Proteome Discoverer v2.1 (Thermo Fisher Scientific) and Mascot server v2.3.02 (Matrix Science) using the SwissProt sequence database for Homo sapiens (www.uniprot.org in November 2016) and the cRAP database (common Repository for Adventitious Proteins, https://www.thegpm.org/crap). Standard filtering of the raw data was conducted as per⁷ to remove contaminants and low abundant PSMs. A maximum of two missing values per TMT experiment was allowed and the data was directly exported from Proteome Discoverer at the PSM level for downstream analysis in R.

The experiments were done in triplicate for each condition and each replicate contains 2 x TMT 10-plex sets (which we denote as “set 1” and “set 2”). Thus, in total we conducted a total of 12 experiments, 6 per condition wherein we had 2 sets per replicate containing a total of 20 fractions/channels per replicate. One TMT channel was removed (TMT tag 126 in replicate 2, set 2 in the unstimulated) due to erroneous labelling of insoluble material during the sample preparation for this specific tag. The “Average Reporter S/N” value was recalculated for the nine remaining channels in the corresponding 10plex set and PSMs with a value less than 9.0 were discarded (please see Data Processing in the Supplementary Information of²).

4.2 PSM level data processing

## Unstimulated data
data("psms_thpLOPIT_unstim_rep1_set1")
data("psms_thpLOPIT_unstim_rep1_set2")

data("psms_thpLOPIT_unstim_rep2_set1")
data("psms_thpLOPIT_unstim_rep2_set2")

data("psms_thpLOPIT_unstim_rep3_set1")
data("psms_thpLOPIT_unstim_rep3_set2")

## 12h post LPS-stimulated
data("psms_thpLOPIT_lps_rep1_set1")
data("psms_thpLOPIT_lps_rep1_set2")

data("psms_thpLOPIT_lps_rep2_set1")
data("psms_thpLOPIT_lps_rep2_set2")

data("psms_thpLOPIT_lps_rep3_set1")
data("psms_thpLOPIT_lps_rep3_set2")

psms <- MSnSetList(
  list(
    psms_thpLOPIT_unstim_rep1_set1,
    psms_thpLOPIT_unstim_rep1_set2,
    psms_thpLOPIT_unstim_rep2_set1,
    psms_thpLOPIT_unstim_rep2_set2,
    psms_thpLOPIT_unstim_rep3_set1,
    psms_thpLOPIT_unstim_rep3_set2,
    psms_thpLOPIT_lps_rep1_set1,
    psms_thpLOPIT_lps_rep1_set2,
    psms_thpLOPIT_lps_rep2_set1,
    psms_thpLOPIT_lps_rep2_set2,
    psms_thpLOPIT_lps_rep3_set1,
    psms_thpLOPIT_lps_rep3_set2
  )
)

## add list names for each dataset
.nams <- paste0("Rep", rep(1:3, each = 2), "_set", rep(1:2, 3))
names(psms) <- c(paste0(.nams, "_unstim"), paste0(.nams, "_lps"))

sapply(psms, nrow)

## Rep1_set1_unstim Rep1_set2_unstim Rep2_set1_unstim Rep2_set2_unstim 
##            44137            41860            39297            37919 
## Rep3_set1_unstim Rep3_set2_unstim    Rep1_set1_lps    Rep1_set2_lps 
##            48199            44802            40994            63844 
##    Rep2_set1_lps    Rep2_set2_lps    Rep3_set1_lps    Rep3_set2_lps 
##            43799            44055            46805            45600

for (i in seq(psms)) {
  par(las = 2, oma = c(10, 1, 1, 1))
  naplot(psms[[i]], col = "black", las = 2, reorderColumns = FALSE, 
         main = names(psms)[i], cex.axis = 2)
}

psms.imputed <- lapply(psms, function(z) impute(z, method = "MinDet"))

psms.imputed.norm <- lapply(psms.imputed, function(z) normalise(z,  "sum"))

prots <- lapply(psms.imputed.norm, function(z) 
  combineFeatures(z, method = "median",
                  groupBy = fData(z)[, "Master.Protein.Accessions"]))

## become combining all protein sets we tidy up the feature labels of the dataset
ll <- paste0(rep(c("unst", "lps"), 1, each = 6), ".r",
             rep(1:3, 1, each = 2), ".s", 1:2)
for (i in seq(prots)) {
  prots@x[[i]] <- updateFvarLabels(prots[[i]], label = ll[i], sep = "_")
}

## combine TMT sets and examine each replicate for each condition
unst_r1 <- filterNA(MSnbase::combine(prots[[1]], prots[[2]]))
unst_r2 <- filterNA(MSnbase::combine(prots[[3]], prots[[4]]))
unst_r3 <- filterNA(MSnbase::combine(prots[[5]], prots[[6]]))
lps_r1 <- filterNA(MSnbase::combine(prots[[7]], prots[[8]]))
lps_r2 <- filterNA(MSnbase::combine(prots[[9]], prots[[10]]))
lps_r3 <- filterNA(MSnbase::combine(prots[[11]], prots[[12]]))


tot_prots <- data.frame("Unstimulated" = c(`Replicate 1` = nrow(unst_r1), 
                                           `Replicate 2` = nrow(unst_r2), 
                                           `Replicate 3`= nrow(unst_r3)),
                        "LPS" = c(`Replicate 1` = nrow(lps_r1), 
                                                 `Replicate 2` = nrow(lps_r2), 
                                                 `Replicate 3` = nrow(lps_r3)))

library("knitr")
library("kableExtra")
kable(tot_prots, caption = "Number of quantified proteins per replicate per condition") %>%
  kable_minimal()

4.3 Protein level data processing

As shown in the above section we quantify between ~4800-5800 proteins per replicate. The below Venn diagrams show the overlap between replicates within conditions. We find 3882 proteins common in the unstimulated hyperLOPIT experiment and 4067 in the 12h-LPS stimulated hyperLOPIT experiment.

## Load R packages required for genetrating Venn diagrams
library("VennDiagram")
library("scales")
library("ggplot2")

venn.diagram(
  x = list(featureNames(unst_r1), featureNames(unst_r2), featureNames(unst_r3)),
  category.names = c("Replicate 1" , "Replicate 2 " , "Replicate 3"),
  filename = "venn_unst_replicates.png",
  main = "Unstimulated hyperLOPIT",
  col=c("#440154ff", '#21908dff', '#fde725ff'),
  fill = c(alpha("#440154ff",0.3), alpha('#21908dff',0.3), alpha('#fde725ff',0.3)),
  fontfamily = "sans",
  cat.fontfamily = "sans",
  main.fontfamily = "sans",
  cat.col = c("#440154ff", '#21908dff', '#fde725ff'), 
  output = FALSE
)

venn.diagram(
  x = list(featureNames(lps_r1), featureNames(lps_r2), featureNames(lps_r3)),
  category.names = c("Replicate 1" , "Replicate 2 " , "Replicate 3"),
  filename = "venn_lps_replicates.png",
  main = "LPS stimulated hyperLOPIT 12h-LPS",
  col=c("#440154ff", '#21908dff', '#fde725ff'),
  fill = c(alpha("#440154ff",0.3), alpha('#21908dff',0.3), alpha('#fde725ff',0.3)),
  fontfamily = "sans",
  cat.fontfamily = "sans",
  main.fontfamily = "sans",
  cat.col = c("#440154ff", '#21908dff', '#fde725ff')
)

Figure 4.1: Overlap of protein groups identified between replicated experiments

(Venn diagrams as per supplementary figure 1 of [1])

4.4 Dataset concatenation

We use the combine function in MSnbase to perform data fusion (concatenation) within each condition to maximise subcellular resolution and thus the reliability in protein classification [as per the pRoloc pipeline.⁶ Trotter et al. ¹⁰ have shown a significant improvement in classifier accuracy when combining replicated experiments in spatial proteomics. Specifically, we concatenate all gradient fractions across replicated experiments to give better discrimination between subcellular niches (as shown in⁶) which gives users the opportunity to resolve niches that may not have been discovered when examining replicates alone.

We note (as explained in⁸) that direct comparisons of individual TMT channels in replicated experiments do not provide an adequate, goal-driven assessment of different experiments due to the nature of the gradient fraction collection, where quantitative channels do not correspond to identical selected fractions along the gradient.

Concatenate replicates for each condition,

## combine replicates for each condition
unst_full <- filterNA(MSnbase::combine(unst_r1, unst_r2))
unst_full <- filterNA(MSnbase::combine(unst_full, unst_r3)) 

lps_full <- filterNA(MSnbase::combine(lps_r1, lps_r2))
lps_full <- filterNA(MSnbase::combine(lps_full, lps_r3)) 

venn.diagram(
  x = list(featureNames(unst_full), featureNames(lps_full)),
  category.names = c("Unstimulated" , "LPS-stimulated"),
  filename = "venn_hyperLOPIT.png",
  main = "Subset of common proteins between hyperLOPIT conditions",
  col=c("#440154ff", '#21908dff'),
  fill = c(alpha("#440154ff",0.3), alpha('#21908dff',0.3)),
  fontfamily = "sans",
  cat.fontfamily = "sans",
  main.fontfamily = "sans",
  cat.col = c("#440154ff", '#21908dff'),
  margin = 0.05,
  cat.pos = c(-140, 140),
  cat.dist = c(0.05, 0.05)
)

Figure 4. Overlap of protein groups identified in both conditions

We find 3288 proteins common all 12 experiments. In the next code chunk we subset our data to keep only proteins common between all experiments.

## keep only those in common between conditions
cmn <- intersect(featureNames(unst_full), featureNames(lps_full))
lps <- lps_full[cmn, ]
unst <- unst_full[cmn, ]

lps <- addMarkers(lps, markers = "csv/markers_2019.28.01.csv")

## Markers in data: 783 out of 3288

## organelleMarkers
##      40S/60S Ribosome             Chromatin               Cytosol 
##                    75                    73                    52 
## Endoplasmic Reticulum       Golgi Apparatus              Lysosome 
##                    69                    39                    49 
##          Mitochondria             Nucleolus               Nucleus 
##                   197                    38                    91 
##            Peroxisome       Plasma Membrane               unknown 
##                    29                    71                  2505

unst <- addMarkers(unst, markers = "csv/markers_2019.28.01.csv")

## Markers in data: 783 out of 3288

## organelleMarkers
##      40S/60S Ribosome             Chromatin               Cytosol 
##                    75                    73                    52 
## Endoplasmic Reticulum       Golgi Apparatus              Lysosome 
##                    69                    39                    49 
##          Mitochondria             Nucleolus               Nucleus 
##                   197                    38                    91 
##            Peroxisome       Plasma Membrane               unknown 
##                    29                    71                  2505

source("r/prettymap.R")
## set organelle colours
mycol <- c("#88CCEE", "#332288", "#53CAB7", "#0170b4", "#204f20", "#990000",
           "#E69F00", "#DDCC77", "#E18493", "#AA4499", "#D55E00")
setStockcol(mycol)
setUnknownpch(16)


## Generate the t-SNE coords (Figure 4 of [1])
set.seed(399)
tsne_unst <- plot2D(unst, method = "t-SNE", plot = FALSE)
set.seed(399)
tsne_lps <- plot2D(lps, method = "t-SNE", plot = FALSE)


## plot the data as t-SNE maps
par(mfrow = c(1, 2))        
oo = c(1:3, 10, 5, 4, 6:7, 9, 8, 11)
myleg <- c(getMarkerClasses(unst), "Unknown")

prettyTSNE(tsne_unst, unst, fcol = "markers", 
          main = "Subcellular markers\nunstimulated THP1-cells",
          orgOrder = oo, mainCol = mycol, oulineCol = darken(mycol))
prettyTSNE(tsne_lps, lps, 
          main = "Subcellular markers\n12h-LPS stimulated THP1-cells",
          orgOrder = oo, mainCol = mycol, oulineCol = darken(mycol))

## add a legend
plot(NULL, xaxt='n',yaxt='n',bty='n',ylab='',xlab='', xlim=0:1, ylim=0:1)
legend("topleft", legend = myleg, col = mycol, bty = "n",
       pch = 19, cex = 1.2, pt.cex = 1.6)

4.5 Classification using a Bayesian framework

Proteins are assigned to one of the 11 subcellular marker classes using TAGM-MCMC.¹² TAGM-MCMC is a semi-supervised Bayesian generative classifier based on a T-Augmented Gaussian Mixture model that uses Bayesian computation performed using Markov-chain Monte Carlo. It is one of the many supervised machine learning methods available in the pRoloc package, but currently the only Bayesian method that is available in the package for the prediction of subcellular location for spatial proteomics data.

The advantage of using TAGM, and Bayesian methods in general, over classic supervised machine learning approaches is the ability for the algorithm to quantify the localisation uncertainty. The TAGM classifier models each annotated subcellular class using a Gaussian distribution and thus the full dataset can be modelled as a mixture of Gaussians. As described in¹² an outlier component probability is also included in this model which is mathematically described as a multivariate student’s t-distribution which gives extra information regarding the

Full details of the TAGM Bayesian methodology and mathematics are found in¹² and a detailed Bioconductor workflow for Bayesian analysis of spatial proteomoics data can be found in.¹³ We follow this workflow for the analysis used here.

## Load the coda library for calculating the Gelman  diagnostics
library("coda")

## Train TAGM-MCMC
p_lps <- tagmMcmcTrain(object = lps,
                       numIter = 25000,
                       burnin = 5000,
                       thin = 10,
                       numChains = 9)

p_unst <- tagmMcmcTrain(object = unst,
                        numIter = 25000,
                        burnin = 5000,
                        thin = 10,
                        numChains = 9)

## All chains look good and all oscillate around an average of 490 outliers
out_unst <- mcmc_get_outliers(p_unst)
out_lps <- mcmc_get_outliers(p_lps)

## Calculate Gelman's Diagnostic
gelman.diag(out_unst)   
gelman.diag(out_lps)

# check all pairs as per the f1000 vignette [10]
gelman.diag(out_unst[1:2])
gelman.diag(out_unst[1:3]) # etc...

mcmc_unst_pooled <- mcmc_pool_chains(mcmc_unst)
mcmc_lps_pooled <- mcmc_pool_chains(mcmc_lps)

mcmc_unst_pooled <- tagmMcmcProcess(mcmc_unst_pooled)
mcmc_lps_pooled <- tagmMcmcProcess(mcmc_lps_pooled)

summary(mcmc_unst_pooled@summary@posteriorEstimates)
summary(mcmc_unst_pooled@summary@tagm.joint)

summary(mcmc_lps_pooled@summary@posteriorEstimates)
summary(mcmc_lps_pooled@summary@tagm.joint)

mcmc_lps <- mcmc_lps_pooled
mcmc_unst <- mcmc_unst_pooled

unst <- tagmMcmcPredict(unst, params = mcmc_unst_pooled, probJoint = TRUE)
lps <- tagmMcmcPredict(lps, params = mcmc_lps_pooled, probJoint = TRUE)

## Load the data from pRolocdata
data("thpLOPIT_unstimulated_mulvey2021")
data("thpLOPIT_lps_mulvey2021")

## We subset for proteins common in all experiments
cmn <- intersect(featureNames(thpLOPIT_unstimulated_mulvey2021),
                 featureNames(thpLOPIT_lps_mulvey2021))
unst <- thpLOPIT_unstimulated_mulvey2021[cmn, ]
lps <- thpLOPIT_lps_mulvey2021[cmn, ]


## Generate the t-SNE coords 
set.seed(399)
tsne_unst <- plot2D(unst, method = "t-SNE", plot = FALSE)
set.seed(399)
tsne_lps <- plot2D(lps, method = "t-SNE", plot = FALSE)

## plot TAGM assignments
par(mfrow = c(1, 2))
mycol <- paste0(getStockcol())
oo <- c(1:3, 10, 5, 4, 6:7, 9, 8, 11)
myleg <- c(getMarkerClasses(unst), "Unknown")

prettyTSNE(tsne_unst, unst, fcol = "tagm.mcmc.allocation", 
          main = "TAGM-MCMC allocation\nunstimulated THP1-cells",
          orgOrder = oo, mainCol = mycol, oulineCol = darken(mycol))
prettyTSNE(tsne_lps, lps, fcol = "tagm.mcmc.allocation",
          main = "TAGM-MCMC allocation\n12h-LPS stimulated THP1-cells",
          orgOrder = oo, mainCol = mycol, oulineCol = darken(mycol))

plot(NULL, xaxt='n',yaxt='n',bty='n',ylab='',xlab='', xlim=0:1, ylim=0:1)
legend("topleft", legend = myleg, col = mycol, bty = "n",
       pch = 19, cex = 1.2, pt.cex = 1.6)

   localisation probability = posterior probability * outlier probability  

fData(unst)$localisation.prob <- fData(unst)$tagm.mcmc.probability * (1 - fData(unst)$tagm.mcmc.outlier)
fData(lps)$localisation.prob <- fData(lps)$tagm.mcmc.probability * (1 - fData(lps)$tagm.mcmc.outlier)

## localisation probability = posterior probability * outlier probability
fData(unst)$localisation.prob <- fData(unst)$tagm.mcmc.probability * (1 - fData(unst)$tagm.mcmc.outlier)
fData(lps)$localisation.prob <- fData(lps)$tagm.mcmc.probability * (1 - fData(lps)$tagm.mcmc.outlier)

## Threshold = (posterior > 0.99) * 1 - (outlier < 1e-6)
t_loc <- 0.99 * (1 - 1e-6)
unst <- getPredictions(unst, 
                       fcol = "tagm.mcmc.allocation", 
                       scol = "localisation.prob", 
                       t = t_loc)

## ans
##      40S/60S Ribosome             Chromatin               Cytosol 
##                    83                   153                   574 
## Endoplasmic Reticulum       Golgi Apparatus              Lysosome 
##                   352                    47                   261 
##          Mitochondria             Nucleolus               Nucleus 
##                   477                    39                   315 
##            Peroxisome       Plasma Membrane               unknown 
##                    34                   165                   788

lps <- getPredictions(lps, 
                      fcol = "tagm.mcmc.allocation", 
                      scol = "localisation.prob", 
                      t = t_loc)

## ans
##      40S/60S Ribosome             Chromatin               Cytosol 
##                    91                   103                   440 
## Endoplasmic Reticulum       Golgi Apparatus              Lysosome 
##                   387                    43                   266 
##          Mitochondria             Nucleolus               Nucleus 
##                   475                    39                   385 
##            Peroxisome       Plasma Membrane               unknown 
##                    40                   227                   792

## plot TAGM assignments
par(mfrow = c(1, 2))
mycol <- paste0(getStockcol())
oo = c(1:3, 10, 5, 4, 6:7, 9, 8, 11)
myleg <- c(getMarkerClasses(unst), "Unknown")

prettyTSNE(tsne_unst, unst, fcol = "tagm.mcmc.allocation.pred", 
          main = "TAGM-MCMC final assignment\nunstimulated THP1-cells",
          orgOrder = oo, mainCol = mycol, oulineCol = darken(mycol))

## Warning in plot.window(...): "oulineCol" is not a graphical parameter

## Warning in plot.xy(xy, type, ...): "oulineCol" is not a graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "oulineCol" is not
## a graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "oulineCol" is not
## a graphical parameter

## Warning in box(...): "oulineCol" is not a graphical parameter

## Warning in title(...): "oulineCol" is not a graphical parameter

prettyTSNE(tsne_lps, lps, fcol = "tagm.mcmc.allocation.pred",
          main = "TAGM-MCMC final assignment\n12h-LPS stimulated THP1-cells",
          orgOrder = oo, mainCol = mycol, oulineCol = darken(mycol))

## Warning in plot.window(...): "oulineCol" is not a graphical parameter

## Warning in plot.xy(xy, type, ...): "oulineCol" is not a graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "oulineCol" is not
## a graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "oulineCol" is not
## a graphical parameter

## Warning in box(...): "oulineCol" is not a graphical parameter

## Warning in title(...): "oulineCol" is not a graphical parameter

plot(NULL, xaxt='n',yaxt='n',bty='n',ylab='',xlab='', xlim=0:1, ylim=0:1)
legend("topleft", legend = myleg, col = mycol, bty = "n",
       pch = 19, cex = 1.2, pt.cex = 1.6)

## Unstimulated protein profiles
par(mfrow = c(4, 3))
cl <- getMarkerClasses(unst)
for (i in seq(cl)) {
  par(mar = c(8, 4, 8, 2), las = 2)
  plotDist(unst[ fData(unst)$markers == getMarkerClasses(unst)[i], ], 
           pcol = paste0(getStockcol()[i], 90), xlab = "", ylim = c(0, .9))
  title(main = cl[i])
}

## LPS protein profiles
par(mfrow = c(4, 3))
cl <- getMarkerClasses(lps)
for (i in seq(cl)) {
  par(mar = c(8, 4, 8, 2), las = 2)
  plotDist(lps[ fData(lps)$markers == getMarkerClasses(lps)[i], ], 
           pcol = paste0(getStockcol()[i], 90), xlab = "", ylim = c(0, .9))
  title(main = cl[i])
}

loc_unst <- getMarkers(unst, "tagm.mcmc.allocation.pred")

## organelleMarkers
##      40S/60S Ribosome             Chromatin               Cytosol 
##                    83                   153                   574 
## Endoplasmic Reticulum       Golgi Apparatus              Lysosome 
##                   352                    47                   261 
##          Mitochondria             Nucleolus               Nucleus 
##                   477                    39                   315 
##            Peroxisome       Plasma Membrane               unknown 
##                    34                   165                   788

loc_lps <- getMarkers(lps, "tagm.mcmc.allocation.pred")

## organelleMarkers
##      40S/60S Ribosome             Chromatin               Cytosol 
##                    91                   103                   440 
## Endoplasmic Reticulum       Golgi Apparatus              Lysosome 
##                   387                    43                   266 
##          Mitochondria             Nucleolus               Nucleus 
##                   475                    39                   385 
##            Peroxisome       Plasma Membrane               unknown 
##                    40                   227                   792

kable(rbind("Unstimulated" = table(loc_unst), "LPS" = table(loc_lps)))

source("R/intersection.R")
## Add a label to the LPS dataset so when we combine it  
## with the unstimulated the information is not lost
.lps <- lps
.lps <- updateFvarLabels(.lps, "LPS")
cmb <- MSnbase::combine(unst, .lps)  # now combine all results

df_all <- compareDatasets(cmb, 
                          fcol1 = "localisation.pred",
                          fcol2 = "localisation.pred.LPS")
ggheatmap(df_all, title = "TAGM locations between conditions")

4.6 Translocations

One of the main aims of the THP-1 study in² was to use the hyperLOPIT technology to investigate the (spatial) changes in protein localisation between the unstimulated and 12h-LPS stimulated conditions. In this section we look at how the assigned location of proteins differ between the two conditions. We examine their assigned location from TAGM and also the difference between their joint posterior probability distributions.

Re-localisations i.e. proteins that do not localise to the same subcellular niche in both conditions, were extracted and categorised as one of four different translocation events:

1. Type 1: organelle-to-organelle - cases where proteins are assigned to different subcellular classes in each condition i.e. those that move from one discrete location to another.
2. Type 2: unknown-to-organelle - cases where proteins in the unstimulated data are labelled as “unknown” but are localised to one of the 11 organelle classes in the 12h-LPS stimulated data.
3. Type 3: organelle-to-unknown - cases where proteins are assigned a location to one of the 11 organelle classes in the unstimulated data but in the 12h-LPS data are labelled as “unknown”.
4. Type 4: unknown-to-unknown - dynamic proteins that were not assigned a subcellular niche. These are proteins which have a difference in their joint probability distribution equal to 1 (as calculated by the natural L2 distance, also known as the Euclidean norm).

The L2 distance is denoted by \[ d_{_{L2}norm}(x, y) = \sqrt {\sum_{i = 1}^{n}{(x_i - y_i)^2}} \] where x and y are the posterior probabilities for for the unstimulated and 12h-LPS stimulated respectively, for each ith class.

As already mentioned, proteins labelled as “unknown” location, describe proteins that did not meet the classifier threshold to be assigned to one of the 11 discrete organelle categories, as defined in our marker set. We not only wish to examine proteins that move between discrete locations (which we call type 1 translocations above) but we also wish to examine cases where in one condition the protein has been given a discrete location, and in another condition it has been described as “unknown”, and therefore not belonging to one of the 11 subcellular classes. It is important to note that in terms of machine learning, the term “unknown” is not a class in the training set, it is simply as term we use in our pipeline to describe proteins that we can say as not confidently assigned to one of the 11 classes.

These cases are still of great interest as although we may not know where the protein exactly resides in this snapshot of time we do know it is not classed to the same location in the other condition. Unknown cases encompass a number of scenarios for example, a protein having a mixed population and moonlighting between organelles, or it be that the protein localised to a subcellular niche that does not appear in our training data. Whichever reason what we are interested in is not only a change in localion but the size of the probability change, as this implies a change in localisation regardless of organelle residency.

As an extra measure of stringency we enforce that in cases 2-4 any protein considered as “unknown” must also have a high outlier probability (as per the below code chunk out = 1e-3) and thus is clearly not a member of any of the 11 classes.

For every protein the natural L2 distance (Euclidean norm) was calculated between the TAGM joint posterior probability distributions and translocation events were further filtered and ranked based on this distance to give a succinct subset of translocations for further analysis and data mining.

In the below code chunk we call the function getTranslocations which extracts the localisations of every protein in each condition then assigns a translocation type as described above and then calculates the L2 distance between the posterior probability distributions. This information is appended to the fData slot of our MSnSet along with all the other machine learning results.

source("R/getmovers.R")
extract_movers <- getTranslocations(unst, lps, 
                                    fcol = "tagm.mcmc.allocation.pred", 
                                    out = 1e-3)
unst <- extract_movers[[1]]
lps <- extract_movers[[2]]

getMarkers(unst, "translocations")

## organelleMarkers
## type1 type2 type3 type4 
##   112    62    49    30

In summary classed by translocation events we find - - 112 type 1 - 62 type 2 - 49 type 3 - 30 type 4 proteins that move following 12 hour stimulation with LPS.

In the next section we can further examine and plot these events using alluvial and circos plots.

4.6.1 Visualisation of translocations

The following code chunks use functions from the circlize and ggalluvial packages (full source code in this repo and at the end of this vignette) to generate circos and alluvial plots respectively. These plots summarise the flow of proteins between organelle locations and the two conditions.

4.6.1.1 Circos plot of protein translocations across organelles

source("R/circos.R")

## generate my colour scheme
circos_cols <- c(getStockcol(), "grey")
orgs <- c(union(getMarkerClasses(unst), getMarkerClasses(lps)), "unknown")
(colscheme <- setNames(circos_cols, orgs))  # check levels consistent 

##      40S/60S Ribosome             Chromatin               Cytosol 
##             "#88CCEE"             "#332288"             "#53CAB7" 
## Endoplasmic Reticulum       Golgi Apparatus              Lysosome 
##             "#0170b4"             "#204f20"             "#990000" 
##          Mitochondria             Nucleolus               Nucleus 
##             "#E69F00"             "#DDCC77"             "#E18493" 
##            Peroxisome       Plasma Membrane               unknown 
##             "#AA4499"             "#D55E00"                "grey"

## get indices of translocating proteins
tl <- !is.na(fData(unst)$translocations)

## set circos plot parameters
par(mar = c(2, 2, 2, 2), cex = .5)
circos.clear()
circos.par(gap.degree = 4)

## plot circos
customChord(unst[tl, ], lps[tl, ], 
            cols = colscheme,  
            diffHeight  = -0.02, 
            transparency = 0.3, 
            link.sort = FALSE,
            labels = TRUE)

Note: labels were optimised in Inkscape

4.6.1.2 Alluvial (river) plot of protein translocations across organelles

source("r/riverplot.R")
library(ggalluvial)
thp_alluvial <- riverplot(unst[tl, ], lps[tl, ], 
                          cols = colscheme,
                          labels = TRUE)
thp_alluvial + ggtitle("Translocations following 12h-LPS stimulation in THP-1 cells") 

Note: labels were optimised in Inkscape

Figure 12. Alluvial (river) plot showing the 253 protein translocations (Fig. 4C of Mulvey et al. 2021

Alluvial (river) plot showing the 253 protein translocations (Figure 4C of²) proteins which were found to move between organelles at 12h-LPS. The number of proteins that are found to translocate to and from each subcellular compartment is denoted next to the labelled alluvial blocks

4.6.2 Spatial map of translocating proteins

The t-SNE maps below also show the visual localisation of translocating proteins.

source("r/foi.R")
par(mfrow = c(1, 2))

## Unstimulated
prettyTSNE_overlay(tsne_unst, unst, 
                   fcol = "localisation.pred", 
                   main = "Type 1, 2, 3, or 4 translocation events", orgOrder =  oo, 
                   mainCol = paste0(lighten(mycol), 30),
                   oulineCol = paste0(lighten(mycol), 30))
highlightOnPlot(object = tsne_unst, pch = 21, cex = 1.7,
                col = mycol[2],  bg = lighten(lighten("grey")),
                foi = type1, lwd = 2,
                args = list(unst))
highlightOnPlot(object = tsne_unst, pch = 24, cex = 1.7,
                col = darken(mycol[7]), bg = lighten(lighten("#f4f4c8")), 
                foi = type2, lwd = 2, 
                args = list(unst))
highlightOnPlot(object = tsne_unst, pch = 21, cex = 1.7,
                col = darken("purple"), bg = lighten(lighten("#f2a6e3")), 
                foi = type3, lwd = 2,
                args = list(unst))
highlightOnPlot(object = tsne_unst, pch = 24, cex = 1.7,
                col = "#1c1c78", bg = lighten(lighten("#5fa8cc")), 
                foi = type4, lwd = 2,
                args = list(unst))

prettyTSNE_overlay(tsne_lps, lps, fcol = "localisation.pred", 
                   main = "", orgOrder =  oo,
                   mainCol = paste0(lighten(mycol), 30), 
                   oulineCol = paste0(lighten(mycol), 30))
## 12h-LPS 
highlightOnPlot(object = tsne_lps, pch = 21, cex = 1.7,
                col = mycol[2],  bg = lighten(lighten("grey")),
                foi = type1, lwd = 2,
                args = list(lps))
highlightOnPlot(object = tsne_lps, pch = 24, cex = 1.7,
                col = darken(mycol[7]), bg = lighten(lighten("#f4f4c8")), 
                foi = type2, lwd = 2, 
                args = list(lps))
highlightOnPlot(object = tsne_lps, pch = 21, cex = 1.7,
                col = darken("purple"), bg = lighten(lighten("#f2a6e3")), 
                foi = type3, lwd = 2,
                args = list(lps))
highlightOnPlot(object = tsne_lps, pch = 24, cex = 1.7,
                col = "#1c1c78", bg = lighten(lighten("#5fa8cc")), 
                foi = type4, lwd = 2,
                args = list(lps))

plot(NULL, xaxt='n',yaxt='n',bty='n',ylab='',xlab='', xlim=0:1, ylim=0:1)
legend("topleft", legend = c("Type 1: organelle-to-organelle", 
                             "Type 2: unknown-to-organelle", 
                             "Type 3: organelle-to-unknown", 
                             "Type 4: unknown-to-unknown"), 
       col = c(mycol[2], darken(mycol[7]), darken("purple"), "#1c1c78"), 
       pt.bg = c(lighten(lighten("grey")), lighten(lighten("#f4f4c8")), 
                 lighten(lighten("#f2a6e3")), lighten(lighten("#5fa8cc"))),
       bty = "n",
       pch = c(21, 24, 21, 24), pt.cex = 2.2, cex = 1.4)

Figure 13. T-SNE dimensionality reduction showing the 253 relocalising proteins from the hyperLOPIT experiments and plotted by Type 1,2,3, or 4 relocalisaton events, which are represented by grey circles, yellow triangles, pink circles and blue triangles, respectively.

df <- makedf(unst[tl, ], lps[tl, ],
            fcol1 = "localisation.pred",
            fcol2 = "localisation.pred",
            mrkCol1 = "markers",
            mrkCol2 = "markers")

## `summarise()` has grouped output by 'x'. You can override using the `.groups` argument.

df <- df[!df$`Number of translocations`==0, ]
names(df) <- c("Unstimulated", "LPS", "Number of translocations")
kable(df)

4.7 A scaffold for proteins of interest

Mulvey et al (2021) show that hyperLOPIT provides a means to robustly classify unannotated proteins to distinct organelles and to investigate protein relocalisation events. The following spatial maps as generated in² highlight interesting findings from the THP-1 hyperLOPIT dataset.

ADD MINI CIRCOS OR RIVER PLOTS FOR ALL OF THESE SPATIAL MAPS?